
program define dyadic
	version 9.0
	if replay() {
		if ("`e(cmd)'" ~= "dyadic") error 301
		dyadic_replay `0'
	}
	else _dyadic `0'
end

program define _dyadic, eclass
	syntax varlist(min=2) [if] [in], Id(varlist min=2 max=2) [Group(varlist) SYMmetric OLS LOGit NOConstant INComplete]

	display _newline

	gettoken y x : varlist

	gettoken igroup1 igroup2 : id

	opts_exclusive "`ols' `logit'"
	if missing("`logit'") local ols "ols"
	if missing("`noconstant'") {
		tempvar cons
		generate byte `cons' = 1
		local x "`x' `cons'"
	}

	marksample touse
	markout `touse' `igroup1' `igroup2' `group' `y' `x', strok

	_rmdcoll `y' `x' if `touse', noconstant
	local x "`r(varlist)'"

	tempname coef var xepex xpx xpxinv glsm
	tempname A1 A2 xacc1 xacc3 xx xxx AA cov_net AACUM XXCUM
	tempvar ehat phat idg idgroup
	tempfile jjj

quietly{
	preserve
	keep if `touse'

	// prepare the data
	count if (`igroup1' == `igroup2')
	// drop diagonal observation (links between the same node)
	if r(N) {
		noisily display as text "note: `r(N)' diagonal observations were not used in calculations"
		drop if (`igroup1' == `igroup2')
	}
	// add the other side of the link when this is not already present
	if ~missing("`symmetric'") {
		keep `touse' `igroup1' `igroup2' `group' `y' `x'
		save `jjj'
		tempname tvar
		rename `igroup1' `tvar'
		rename `igroup2' `igroup1'
		rename `tvar' `igroup2'
		append using `jjj'
	}
	// "group" is a variable that identifies separate networks: if missing the whole db is a single network
	if missing("`group'") {
		tempvar group
		generate byte `group' = 1
	}
	// LOGIT Estimation
   	if missing("`ols'") {
		logit `y' `x', noconstant robust
		predict `phat'                 		// Phat is exp(XB)/(1+exp(XB))
		generate double `ehat' = `y' - `phat'   // The residual
	}
	// OLS Estimation
	else {
		regress `y' `x', noconstant robust
		predict double `ehat', residuals
	}
	
	matrix `coef' = e(b)
	matrix `var'  = e(V)

  	local nobs = e(N)
	local K    = e(df_m)
 	local invN = 1/`nobs'

	// address the problem of non-complete networks 
	if ~missing("`incomplete'") {
		tabulate `group'
		local totgroup=r(r)
		egen `idg' = group(`group')
		qui: count
		local before   `r(N)'
		local complete    0
		forvalues r = 1/`totgroup' {
			qui: count if `idg' == `r'
			local ngroup = (1+sqrt(1+4*`r(N)'))/2
						
			if !mod(`ngroup',1) {		
				local complete = `complete' + 1
				noisily display as text _n "Note: network " `r' " is complete: option INComplete does not affect it."
			}
		}
		if `complete' == `totgroup' {
				noisily display as text _n "Note: all networks are complete: option INComplete ignored."
		}
		else {
			tempfile backup id newid
			tempvar  match n 

			save `backup'
			
			sort `igroup1' `group' 
			collapse (count) `n' = `igroup2', by(`igroup1' `group')
			sort  `igroup1' `group'
			save `id'

			use `backup', clear
			sort `igroup2' `group' 
			collapse (count) `n' = `igroup1', by(`igroup2' `group')

			rename `igroup2' `igroup1'
			merge 1:1 `igroup1' `group' using `id', gen(`match')

			sum `match'

			if `r(min)' < 3 {                       
				forvalues net = 1/`totgroup' {	
					count if `match' == 2 & `group' == `net'
					if `r(N)' == 1 noisily display as text _n "There is "  `r(N)' " individual on network "  `net' " that never appears as `igroup2'."
					if `r(N)' >  1 noisily display as text _n "There are " `r(N)' " individuals on network " `net' " that never appear as `igroup2'."
					count if `match' == 1 & `group' == `net'
					if `r(N)' == 1 noisily display as text _n "There is "  `r(N)' " individual on network "  `net' " that never appears as `igroup1'."
					if `r(N)'  > 1 noisily display as text _n "There are " `r(N)' " individuals on network " `net' " that never appear as `igroup1'."
				}
				
				drop if (`match' == 3)
				gen `igroup2' = `igroup1'
				save `newid', replace
				
				use `backup', clear
				append using `newid'
			}			

			else {                                  
				use `backup', clear
			}

			
			keep `group' `igroup1' `igroup2' `y' `x' `ehat'
			reshape wide `y' `x' `ehat', i(`igroup1' `group') j(`igroup2')
			reshape long `y' `x' `ehat', i(`igroup1' `group') j(`igroup2')

			reshape wide `y' `x' `ehat', i(`igroup2' `group') j(`igroup1')
			reshape long `y' `x' `ehat', i(`igroup2' `group') j(`igroup1')

			drop if (`igroup1' == `igroup2')
			recode `y' `x' `ehat' (.=0)
			qui: count
			if `r(N)' - `before' == 1 {
				noisily display as text _n "Note: there is " `r(N)' - `before' " observation missing to complete the network(s)."
			}
			else {
				noisily display as text _n "Note: there are " `r(N)' - `before' " observations missing to complete the network(s)."
			}
		}
	}
	
	foreach i of local x {
		tempvar tmp
		generate double `tmp'=`i'*`ehat'
		local ee "`ee' `tmp'"
	}

	tempvar rank merge

		sort `group' `igroup2' `igroup1'
		generate `rank'=_n
		save `jjj', replace

		foreach i of local ee {
			tempname tmp
			rename `i' `tmp'
			local e "`e' `tmp'"
		}

	sort `group' `igroup1' `igroup2'
	replace `rank'=_n

	merge `rank' using `jjj', sort _merge(`merge') keep(`ee')

	matrix `AACUM'=J(`K',`K',0)
	matrix `XXCUM'=J(`K',`K',0)

	tabulate `group'
	local totgroup=r(r)
	egen `idgroup'=group(`group')

	forvalues ii=1/`totgroup' {
		local c =cond(mod(`ii',50)==0, "50", "_c")
		noisily display as text "." `c'

		count if `idgroup'==`ii'

		// count the number of observations in the group 
		// (make N # individuals and x number of obs, x = N(N-1). Knowing x (=`r(N)') N is correctly calculated with ngroup
		local ngroup=(1+sqrt(1+4*`r(N)'))/2

		if mod(`ngroup',1) {
			noisily display as error _n "network matrix for the group " `ii' " is not square (possibly because of missing observations)." _n "Try using the option INComplete."
			exit 503
		}

		// must adjust this so that only within group
		matrix accum `xepex'=`e' if `idgroup'==`ii', noconstant
		matrix accum `xpx'  =`x' if `idgroup'==`ii', noconstant

		// first the direct correlations
		matrix define `glsm' =J(`ngroup',`ngroup',1)

		sort `igroup1' `igroup2'
		matrix glsaccum `xx'=`e' `ee' if `idgroup'==`ii', group(`igroup1') row(`igroup2') glsmat(`glsm') noconstant

		matrix `A1'=`xx'[1..`K',1..`K']            // this select the top-left K x K matrix in xx
		matrix `A2'=`xx'[`K'+1...,`K'+1...]        // this select the bottom-right K x K matrix in xx
	 noi
		matrix `xacc1'=`xx'[1..`K',`K'+1...]
  		matrix `xacc1'=0.5*(`xacc1'+`xacc1'')

		matrix accum `xxx'=`e' `ee' if `idgroup'==`ii', noconstant
		matrix `xacc3'=`xxx'[1..`K',`K'+1...]

		// now put all the pieces together
		matrix `AA'=`A1'+`A2'-`xepex'+2*`xacc1'-`xacc3'

		matrix `AACUM'=`AACUM'+`invN'*`AA'
		matrix `XXCUM'=`XXCUM'+`xpx'

	} // forvalues

	if missing("`ols'") {
		tempvar sqdphat
		generate double `sqdphat' = sqrt(`phat'*(1-`phat'))

		foreach i of local x {
        	tempvar tmp
        	generate double `tmp'=`i'*`sqdphat'
			local stf "`stf' `tmp'"
        }
		tempname dg
        matrix accum `dg'=`stf', noconstant
		matrix 		 `dg'= `invN'*`dg'

		matrix `cov_net' = `invN'*inv(`dg'*inv(`AACUM')*`dg'')
	}
	else {
		matrix `xpxinv'=`nobs'*syminv(`XXCUM')
		matrix `cov_net'=1/(`nobs'-`K')*`xpxinv'*`AACUM'*`xpxinv'
	}
} // quietly
	restore

	*local x : subinstr local x "`cons'" "_cons"
	matrix rownames `cov_net' = `x'
	matrix colnames `cov_net' = `x'
	matrix rownames `var'     = `x'
	matrix colnames `var'     = `x'
	matrix colnames `coef'    = `x'

	ereturn post `coef' `var', depname("`y'") esample(`touse') obs(`nobs')

	ereturn scalar groups = `totgroup'

	local method =cond(missing("`ols'"), "logistic", "linear")

	ereturn local title "Grouped Dyadic `method' regression"
	ereturn local cmd "dyadic"
	ereturn local model "`ols'`logit'"
	ereturn local vcetype "Robust"
	ereturn local depvar "`y'"
	
	matrix `var'=e(V)
   	ereturn matrix white_var=`var'
	ereturn repost V=`cov_net'

   	dyadic_replay `0'
end

program define dyadic_replay

	display _newline
	display as text e(title) _c
  	display as text _col(40) "Number of obs"    _col(57) "= " as res %7.0f e(N)
  	display as text _col(40) "Number of groups" _col(57) "= " as res %7.0f e(groups)
  	display ""

	 tempname Tab
	.`Tab' = ._tab.new, col(5) lmargin(0) ignore(.b)
	// column        1      2     3     4     5
	.`Tab'.width	13    |12    14    16     9
	.`Tab'.strfmt    .   %11s     .     .     .
	.`Tab'.pad       .      2     2     6     1
	.`Tab'.numfmt    .  %9.0g %9.0g %9.0g %8.2f
	.`Tab'.strcolor  . result    .     .     .

    .`Tab'.sep, top
	.`Tab'.titlefmt  .      .  %10s      .     .
	.`Tab'.titles   "" "" "Robust" "" ""
	.`Tab'.titlefmt  .   %11s  %12s  %16s   %7s
	.`Tab'.titles	"`e(depvar)'" "Coef." "Std. Err." "Gr.Dyadic s.e." "t"
	.`Tab'.sep

	tempname b se se1 t cov_net white_var
*	matrix `cov_net' = e(cov_net)
	matrix `cov_net' = e(V)
	matrix `white_var' = e(white_var)

	foreach i in `: colnames e(b)' {
		scalar `b'=_b[`i']
*		scalar `se'=_se[`i']
	   	matrix `se'=`white_var'["`i'","`i'"]
		scalar `se'=`se'[1,1]
		scalar `se'=sqrt(`se')

	   	matrix `se1'=`cov_net'["`i'","`i'"]
		scalar `se1'=`se1'[1,1]
		scalar `se1'=sqrt(`se1')
		scalar `t'=`b'/`se1'
		local na = abbrev("`i'",12)

		.`Tab'.row "`na'" `b' `se' `se1' `t'
	}
	.`Tab'.sep, bottom
end
